#include "StdAfx.h"
#include "Fourier.h"

#include <complex>
#include <assert.h>

COLORREF Fourier::ColorAtAngle(float inAngle) const
{
	using namespace DirectX;
	XMVECTOR color = ColorComponentsAtAngle(inAngle);
	color = XMVectorClamp(color, XMVectorZero(), XMVectorSet(255, 255, 255, 0));

	int r = ::lround(XMVectorGetX(color));
	int g = ::lround(XMVectorGetY(color));
	int b = ::lround(XMVectorGetZ(color));
	return RGB(r, g, b);
}

DirectX::XMVECTOR Fourier::ColorComponentsAtAngle(float inAngle) const
{
	using namespace DirectX;
	return mA0 + ColorComponentsWithoutConstantAtAngle(inAngle);
}

int Fourier::GetDegree() const
{
	return mDegree;
}

DirectX::XMVECTOR Fourier::GetConstantCoefficient() const
{
	return mA0;
}

DirectX::XMVECTOR Fourier::IntegralWithoutConstantAtAngle(float inAngle) const
{
	using namespace DirectX;
	XMVECTOR result = XMVectorZero();
	auto e1 = exp(std::complex<float>(0, inAngle));
	auto e = e1;
	for (int i = 0; i < mDegree; ++i)
	{
		auto scaledExponent = e / float(i + 1);
		result += mA[i] * scaledExponent.imag() - mB[i] * scaledExponent.real();
		e *= e1;
	}
	return result;
}

DirectX::XMVECTOR Fourier::ColorComponentsWithoutConstantAtAngle(float inAngle) const
{
	using namespace DirectX;
	XMVECTOR result = XMVectorZero();
	auto e1 = exp(std::complex<float>(0, inAngle));
	auto e = e1;
	for (int i = 0; i < mDegree; ++i)
	{
		result += mA[i] * e.real() + mB[i] * e.imag();
		e *= e1;
	}
	return result;
}

DirectX::XMVECTOR Fourier::DerivativeAtAngle(float inAngle) const
{
	using namespace DirectX;
	XMVECTOR result = XMVectorZero();
	auto e1 = exp(std::complex<float>(0, inAngle));
	auto e = e1;
	for (int i = 0; i < mDegree; ++i)
	{
		auto scaledExponent = e * float(i + 1);
		result += mB[i] * scaledExponent.real() - mA[i] * scaledExponent.imag();
		e *= e1;
	}
	return result;
}

float Fourier::SquareOfDerivativeAtAngle(float inAngle) const
{
	using namespace DirectX;
	XMVECTOR derivative = DerivativeAtAngle(inAngle);
	return XMVector3Dot(derivative, derivative).m128_f32[0];
}

float Fourier::IntegralOfSquareAroundCircle() const
{
	float sumOfSquares = 0;
	using namespace DirectX;
	for (int i = 0; i < mDegree; ++i)
	{
		sumOfSquares += XMVector3Dot(mA[i], mA[i]).m128_f32[0] +
			XMVector3Dot(mB[i], mB[i]).m128_f32[0];
	}
	return XM_PI * sumOfSquares;
}

//-------------------------------------------------------
FourierCalculator::FourierCalculator()
	: mDegree(0)
{
}

void FourierCalculator::SetDegree(int inDegree)
{
	if (inDegree == mDegree)
	{
		return;
	}
	assert(inDegree <= kMaxDegree);
	mDegree = inDegree;
	mNumberOfIterations = 2 * mDegree + 1;

	mCosines[0] = 1;
	mSines[0] = 0;
	
	float phi = float(2 * M_PI / mNumberOfIterations);
	auto e1 = exp(std::complex<float>(0, phi));
	auto e = e1;
	for (int i = 1, j = mNumberOfIterations - 1; i < j; ++i, --j)
	{
		float oldC = mCosines[j] = mCosines[i] = e.real();
		mSines[j] = -(mSines[i] = e.imag());
		e *= e1;
	}
}

void FourierCalculator::BeginUpdatingFourier(Fourier& outFourier)
{
	outFourier.mDegree = mDegree;
	mCurrentIterationIndex = -1;

	using namespace DirectX;
	outFourier.mA0 = XMVectorZero();
	for (int i = mDegree; i--;)
	{
		outFourier.mA[i] = outFourier.mB[i] = XMVectorZero();
	}
}

bool FourierCalculator::GetNextPoint(float& outX, float& outY)
{
	int i = ++mCurrentIterationIndex;
	if (i < mNumberOfIterations)
	{
		outX = mCosines[i];
		outY = mSines[i];
		return true;
	}
	return false;
}

void FourierCalculator::AddColorToFourier(DirectX::XMVECTOR inColor,
	Fourier& ioFourier)
{
	using namespace DirectX;
	ioFourier.mA0 += inColor;
	int index = 0;
	for (int i = 0; i < mDegree; ++i)
	{
		index += mCurrentIterationIndex;
		if (index >= mNumberOfIterations)
		{
			index -= mNumberOfIterations;
		}
		ioFourier.mA[i] += inColor * mCosines[index];
		ioFourier.mB[i] += inColor * mSines[index];
	}
}

void FourierCalculator::EndUpdatingFourier(Fourier& outFourier)
{
	using namespace DirectX;
	float coef = 1.0f / mNumberOfIterations;
	outFourier.mA0 *=  coef;

	coef *= 2.0f;
	for (int i = mDegree; i--;)
	{
		outFourier.mA[i] *= coef;
		outFourier.mB[i] *= coef;
	}
}
